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For axially symmetric solutions of Einstein equations there exists a gauge which has the remark- 
able property that the total mass can be written as a conserved, positive definite, integral on the 
spacelike slices. The mass integral provides a nonlinear control of the variables along the whole 
evolution. In this gauge, Einstein equations reduce to a coupled hyperbolic-elliptic system which is 
formally singular at the axis. As a first step in analyzing this system of equations we study linear 
perturbations on flat background. We prove that the linear equations reduce to a very simple sys- 
tem of equations which provide, thought the mass formula, useful insight into the structure of the 
full system. However, the singular behavior of the coefficients at the axis makes the study of this 
linear system difficult from the analytical point of view. In order to understand the behavior of the 
solutions, we study the numerical evolution of them. We provide strong numerical evidence that the 
system is well-posed and that its solutions have the expected behavior. Finally, this linear system 
allows us to formulate a model problem which is physically interesting by itself, since it is connected 
with the linear stability of black holes solutions in axial symmetry. This model can contribute 
significantly to solve the nonlinear problem and at the same time it appears to be tractable. 

PACS numbers: 04.20.Ex, 04.25. D, 02.30.Jr 



I. INTRODUCTION 



Axisymmetric spacetimes has been studied mainly for 
two reasons. The first one is that they often appear in 
astrophysical models like rotating stars and black holes. 
The second is because in the presence of any symmetry 
Einstein equations simplify considerable and hence these 
spacetimes are useful as intermediate step to understand 
more complex problems. In particular, axially symmetric 
gravitational waves in vacuum do not carry angular mo- 
mentum, this represents an important simplification in 
the dynamics. Also, axial symmetry is the only symme- 
try compatible with asymptotic flatness and non-trivial 
gravitational radiation [l]. From this perspective, axially 
symmetric gravitational waves are the simplest possible 
waves emitted from isolated sources. And hence they 
represent the natural candidates to study the strong field 
dynamics of gravitational waves in Einstein equations. 

However, axial symmetry presents a major difficulty. 
To take advantage of the symmetry an adapted coordi- 
nate system should be used in order to reduce the field 
equations to a lower-dimensional system (there is a well 
known procedure to do this for any symmetry in a ge- 
ometrical way [2|, we review this result in Sec. MI A[) . 
The problem is that the norm of the axial Killing vector 
vanishes at the axis, and hence the reduced equations are 
formally singular there. 



This difficulty is so severe that until recently axially 
symmetric spacetimes have not been studied in detail 
even using numerical techniques (see chapter 10.4 in Q 
and references therein). In a number of recent articles 
[3 , 0] j , , @| this kind of singular behavior has been 
successfully implemented numerically. There is however 
no analytical study of axial symmetry in the dynamical 
regime (see the review article Q for results for other kind 
of symmetries). In fact, it can be argued that this sin- 
gular behavior near the axis is so complicated that the 
axially symmetric case is as hard as the full general case 
from the analytical point of view. 

There exists however a new ingredient that makes, in 
our opinion, the problem worth studying. In the article 
[Tfj| it has been proved that there exists a gauge in axial 
symmetry such that the total mass of the spacetime can 
be written as a positive definite volume integral over the 
spacelike slices of the foliation. Moreover, this integral 
is conserved along the evolution. This conserved integral 
control the norm of the fields along the whole evolution. 
This is certainly a very desirable property of this gauge 
which is not present in the general, non-symmetric, case. 
Also, this mass integral formula appears to be connected 
with stability properties of black holes in axial symmetry 

n 

The gauge mentioned above is a combination of 
the well known maximal condition for the lapse and 
the choice of isothermal coordinates (also called quasi 
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isotropical) for the shift. The later condition is only 
possible in axial symmetry. We call it the maximal- 
isothermal gauge. This gauge has been known for long 
time (see [l2[ and [l3j]) but without noticing this prop- 
erty of the mass. It is also important to emphasize that 
this gauge is the one used in most of the recent numer- 
ical computations 0] 0] [H [H (examples of other gauge 
choices in axial symmetry are given in [3 0). That 
is, this gauge has not only desirable analytical properties 
but it is also useful for numerical studies. 

The very basic question of well-posedness of the equa- 
tions in this gauge is open. This question is rather subtle 
because of the singular behavior mentioned above. The 
standard theory in partial differential equations does not 
seems to apply in a direct way. This is the problem we 
want to study in this article. In order to do this, the first 
step is to study the linearization of the equation around 
fixed solutions. We chose Minkowski as a background for 
simplicity. As we describe in the next section, we obtain 
a remarkable simple system of linear equation together 
with a conserved quantity which corresponds to the mass 
of the spacetimc up to second-order corrections. This sys- 
tem allows us to formulate the problem of well-posedness 
in a simplified setting which is nevertheless relevant and 
physically interesting. Remarkably enough, even for this 
linear system the well-posedness appears to be a nontriv- 
ial problem. In order to get insight into this problem we 
numerically evolve these equations to provide evidences 
that the system is in fact well-posed and that the solu- 
tions have the expected behavior. 

If the local existence problem is so complicated in this 
gauge one can wonder what can be said about the global 
behavior of the evolution, which is, of course, the ulti- 
mate goal. However, many of the main complications 
of this gauge are already present in the well-posedness 
problem because they are related with the local behav- 
ior of the fields at the symmetry axis. If one can solve 
them at the linearized level in a satisfactory way there 
is a good chance that the mass integral formula can be 
used to control the global evolution in some way. Also, 
the well-posedness of the linear equations are relevant by 
themselves for the following two reasons. First, the mass 
formula at the linear level can in principle be used to 
prove linear stability in axial symmetry of a background 
solution like a black hole. Second, the well-posedness of 
the linear equations and the mass formula give insight on 
appropriate boundary conditions on a bounded domain. 
In particular, the mass formula allows us to calculate the 
gravitational waves that leave or enter a bounded do- 



their main properties. In particular, in this section we 
discuss the mass conservation and boundary conditions 
on a bounded domain. In Sec. I VII we describe the nu- 
merical techniques used to evolve these equations. And 
in Sec. I VIII we present the numerical results. Finally, in 
Sec. IVIIII we conclude with a discussion of the relevant 
open problems. 



II. MAIN RESULTS 

This article has two main results. The first one is to 
prove that the linearized Einstein vacuum equations in 
the maximal-isothermal gauge reduce to a very simple 
set of equations together with a conserved quantity. This 
conserved quantity is the mass up to second-order correc- 
tions and it is written as a positive definite integral over a 
spacelike surface, which has a similar form as the energy 
of the wave equation. This property of the mass, which 
only holds in this gauge, is of course what distinguished 
this system of equations from any other linearization. 

The second result is the numerical study of these equa- 
tions, together with the analysis of appropriate boundary 
conditions on a finite grid. 

Let us describe the first result. In axial symmetry, 
the dynamical degrees of freedom of the vacuum gravi- 
tational field are prescribed by two functions, which can 
be chosen to be the norm and twist potential of the axial 
Killing vector (see Sec. IIIip . We make, for simplicity, the 
extra assumption that the twist is zero (although we dis- 
cuss the full non-linear equations with twist in Sec. IIII[) . 
This assumption simplify the equations but it is by no 
means essential. In the maximal-isothermal gauge, the 
linearized Einstein equations with respect to a Minkowski 
background reduce to the following two equations for the 
functions v and (3 P (the reason for the notation for the 
last function is that it represents the p component of the 
shift vector as we will see below) 
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These equations are deduced in Sec. IIV1 We have chosen 
cylindrical coordinates (i, p, z). The relevant domain for 
these equations is the half plane p > 0, — oo < z < oo 
denoted by R^_. A dot denotes time derivative and A is 
the flat Laplacian in 2-dimcnsions 



The plan of the article is the following. In Sec. |TT] 
we summarize our mains results. In Sec. IIIII we review 
the axially symmetric, vacuum, Einstein equations. Al- 
though this is well known, the way we process to obtain 
the final equations in the maximal-isothermal gauge is 
slightly different than the standard one used in the nu- 
merical works mentioned above. In Sec. IIVI we derive 
our main linear equations and in Sec. [V] we describe 



Av = dlv 



dlv. 



(3) 



The boundary condition for Eqs. (JTJ) and ^ arise from 
the regularity of the spacetime metric at the axis and the 
standard asymptotically flat fall-off behavior at infinity. 
We discuss this in detail in Sec. IIIIDI and [V] Let us 
present here a summary. Eq. @ is an elliptic equation 
for /3 P , we need to prescribe boundary conditions on . 



3 



On the axis p = we require 

/3 P U=o = 0, (4) 
and at infinity we impose 

pP = Oir- 1 ), (5) 

where r = yj p 2 + z 2 . With these boundary conditions, 
Eq. ([2]) has a unique solution. Eq. ([I]) is a wave equation 
for v, we need to prescribe initial conditions, which are 
functions f(z,p) and g(z,p) such that 

v\ t =o = f, v\t=o = 9- (6) 

The axis represents a timelike boundary for the wave 
equation |T]) , and hence we need to prescribe also bound- 
ary conditions there. This is the delicate part, because 
the equations are singular at the axis and hence we are 
not free to chose arbitrary boundary conditions there. 
From the axial regularity of the spacetime metric we de- 
duce that the initial data / and g should vanish at the 
axis, namely 

/| P =o=0, 5 |p=o = 0. (7) 

In Sec. El using series expansions, we prove that condi- 
tions (0) on the initial data imply that 

v\ p=0 = 0, d p v\ p=0 = 0, (8) 

for all times. Moreover, solutions v and j3 p of equations 
([1]) and |2]) satisfy a parity conditions, namely v is an even 
function of p and (3 P is an odd function of p. These parity 
conditions imply that the spacetime metric is smooth at 
the axis. It is important to emphasize that these con- 
ditions are consequences of equations (Q} and ([2]) alone, 
without any extra requirement. 

In the numerical implementation, Eqs. ([5]) are used as 
boundary conditions at the axis. There are various ways 
to re-express (TTJ and @ in order write conditions ([5]) as 
proper timelike boundary conditions (e.g. Dirichlet or 
Neumann). For example, following [J], in Sec. IVII we 
write them in terms of the rescaled variable v — v /p. 

We are interested in asymptotically flat solutions of (Q]) 
and ([2]). We will argue in Sec. [V] that the typical fall off 
behavior as r — > oo for this kind of solutions is 

v = 0(r- 2 ). (9) 

That is, if we chose initial data / and g which satisfy ([9]) 
then v will satisfy (O for all times. 

All the other components of the linear perturbation 
can be calculated in terms of v and (3 P as follows. In 
our gauge the four dimensional coordinates are given by 
(t,p,z,(f>). A general twist free linear perturbation is 
written as follows 

7 = (cr + 2q)(dp 2 + dz 2 ) + 2ji p dpdt + 2f3 z dzdt + p 2 ad(/) 2 . 

(10) 



Where the functions cr, q, /3 P and f3 z depends only on 
(t, p, z). The function (3 P is given by the other func- 
tions are calculated in terms of v as follows. The func- 
tions q is a time derivative of v 

q = v. (11) 

The function a is determined by the following elliptic 
equation 

( 3 >A<7 = -At>, (12) 
where ^ 3 ^A is defined as 

( 3 )Aa = Aa+^. (13) 
P 

This operator, which appears frequently in the rest of 
the article, is the flat Laplace operator in 3-dimcnsions 
written in cylindrical coordinates and acting on axially 
symmetric functions. The boundary condition for equa- 
tion (fT2]) at the axis is given by 

d p a\ p=0 = 0, (14) 

and at infinity we impose 

(i = 0{r- 1 ). (15) 

Eq. (jT2|) can be also viewed as an equation in R 3 . In this 
case we do not need to prescribe any boundary condition 
at the axis. Condition (f!4| will be automatically satisfied 
for any regular solution. 

Finally, the other component of the shift vector is de- 
termined by the following equation 

A/3* = -2\, (16) 
P 

with boundary condition at the axis 

0„j8"U=o = O, (17) 
and decay condition at infinity 

Z = 0{r- r ). (18) 

The total mass of the system is given by the following 
integral 

m = YeL ( 4 ^~ + {Av)2 + |9ff|2 ) pdpdz ' (19) 

Note that in order to compute the mass we need the 
function cr, which satisfies Eq. (|12p . This equation is 
uncoupled with equations (TTJ and @. The integral ([T9")) 
is conserved. That is, for every solution of |T]) and @ 
which satisfy the boundary conditions (j4j) , ([5]) , (J8|) and 
decay at infinity like (J9j) we have 

m = 0. (20) 
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The conservation law (|20[) is deduced from a local conser- 
vation formula which involves the integrand of the mass 
formula (fT9| . This local conservation law can be also 
used to compute the gravitational waves entering or leav- 
ing a bounded domain. We discuss this in Sec. [V] 

The second main result of this article is the numeri- 
cal study of the system ([T]) and @ . We describe this in 
detail in Sec. IVHl Let us briefly summarize these re- 
sults. The system (TTJ and © appears, from the numer- 
ical evidences, to be well posed and numerically stable. 
In particular, this imply that the functions v and (3 P re- 
main bounded for all times by a constant that depends 
only on the initial data. This is consistent with the linear 
stability of Minkowski spacetime. 

The numerical calculations are, of course, performed 
on a finite grid. Hence, we need to prescribe bound- 
ary conditions on a bounded domain. These conditions 
should be compatible with asymptotically flatness in the 
following sense. Assume we have a sequence of bounded 
domains such that in the limit they cover the half plane 
R^_. If we solved the equations for this sequence of do- 
mains we should recover in the limit the asymptotically 
flat solution described above. There exists many different 
boundary conditions that have this property. In particu- 
lar, homogeneous Dirichlet conditions for (3 P and v. For 
each bounded domain the mass is not conserved. How- 
ever, as the size of the domain increase we expect that 
the mass approach a time independent constant. This is 
precisely what we observe in our numerical calculations. 

For our present goal, this kind of asymptotically flat 
boundary conditions is all what we need. There is, how- 
ever, an interesting extra point here. To model an iso- 
lated system on a finite grid it is important to prescribe 
boundary conditions such that the gravitational radia- 
tion leaves the domain. In general, this is a very dif- 
ficult problem since it is not even clear what we mean 
by gravitational radiation at a finite distance. However, 
as we mention above, in our gauge the mass formula al- 
low us to compute gravitational radiation on a bounded 
domain. Although it appears not to be possible to pre- 
scribe boundary conditions such that the gravitational 
waves always leave the domain, the mass formula sug- 
gests a particular kind of boundary conditions that has 
this behavior in our numerical calculations. That is, un- 
der these boundary conditions, the mass on a bounded 
domain is monotonically decreasing with time for the par- 
ticular kind of initial data used in the computations. We 
emphasize however that we have not been able to prove 
this analytically. We explore this in detail in Sec. IVland 

Em 



III. AXISYMMETRIC VACUUM EINSTEIN 
EQUATIONS 

The purpose of this section is to write the vacuum Ein- 
stein equations for axially symmetric spacetimes in the 
maximal-isothermal gauge. This involves three clearly 



distinguished steps. In the first one, described in Sec. 
IIII Al we perform a symmetry reduction of Einstein equa- 
tions to obtain a set of geometrical equations in the 3- 
dimensional quotient manifold. These equations can be 
viewed as 3-dimensional Einstein equations coupled with 
effective matter sources. In the second step (Sec. IIII B|) 
we chose an arbitrary spacelike foliation in the quotient 
manifold and split the equations in time plus space. In 
Sec. IIII CI we fixes the foliation and the coordinate sys- 
tem. We also write the mass formula in this gauge. Fi- 
nally, in Sec. IIIIDI we discuss boundary conditions at the 
axis and at infinity. 

A. Symmetry reduction 

In this section we perform the symmetry reduction of 
the field equations. We follow @ and See also [f| 

Consider a vacuum solution of Einstein's equations, 
i.e., a four dimensional manifold M with metric g^ (with 
signature ( — h ++)) such that the corresponding Ricci 
tensor vanishes 

(4) 7V = 0. (21) 

Suppose, in addition, that the metric g^ v admits a Killing 
field rf , that is rf satisfies the equation 

V (jU ^) = 0, (22) 

where V M is the connection with respect to g^ v . Greek 
indices /i, ^, • • • denote four dimensional indices. 

We define the square of the norm and the twist of rj^, 
respectively, by 

V = r^rfg^, = e^if^V- (23) 
Using the field Eq. (f2"Tj) it is possible to prove that 

V\„u v] = 0, (24) 

and hence w M is locally the gradient of a scalar field lu 

lu^ = \7fj_uj. (25) 

Let Af denote the collection of all trajectories of rj^, 
and assume that it is a differential 3-manifold. We define 
the metric h^ v on M by 

V9ix.v = h fiv +ri fl ri v . (26) 

The vacuum field equations (|21[) can be written in the 
following form on J\f 

Urj = -(y a v V a r, - V a coV a io), (27) 
V 

Uoj = -V a u;V a r), (28) 
V 

{3) Rab = ^(Va^Vftr? + V a wV fc w). (29) 
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where V a and ( 3 'R a b are the connexion and the Ricci 
tensor of h ab , we have defined □ = V a V a and Latin 
indices a,b . . . denote three dimensional indices on Af. 

Note that the definition of the metric (J^SJ) involves a 
conformal rescaling with respect to the canonical metric 
h Vil defined by 



g vl i =h Vfl + ri 1 rj l/ r] lx . 



That is, we have 



L flV 



Vhuv 



(30) 



(31) 



This rescaling simplify considerably the field equations. 
In particular, on the right hand side of Eq. 1|29[) there 
are no second derivatives of the fields rj and lo (compare, 
for example, with equation (20) in [l(|). 

Finally, we note that Eq. (f27j) can be written in the 
following form 



□E = 
where we have defined 



2 ' 



E = log 77. 



(32) 



(33) 



Up to this point, the only assumption we have made 
is that the spacetime admits a Killing vector field if" 
and that rf is not null, otherwise the metric h ab is not 
defined. If the Killing field is timelike (r) < 0) then the 
metric h a b is Riemannian and the equations ([27|) - ([29]) are 
the stationary Einstein vacuum equations. On the other 
hand, when the Killing vector is spacelike (rj > 0), the 
metric h ab is a is a 3-dimensional Lorenzian metric (we 
chose the signature ( — h +)). In axially symmetry, the 
Killing vector rf 1 is spacelike and its norm vanishes at 
the axis of symmetry. Hence, the equations are formally 
singular at the axis. This singular behavior at the axis 
represents the main difficulty to handle these equations. 

In the Lorenzian case, Eq. (|2"9")l has the form of Ein- 
stein equations in three dimensions, with effective matter 
sources produces by r\ and lo. The effective matter Eqs. 
(f2"T)) - (f2"5|) imply that the energy-momentum tensor de- 
fined in terms of rj and lo by 

1 



Tab = — (V a 7?V677 + VaWV{,w) 

1 



2?7 2 



4// 



2 /iab(V c T ? V c 77 + V c wV c w) , (34) 



is divergence free, i.e. V a T ab = 0. 

A particularly relevant special case is when lo = 0. In 
that case Eqs. (|27 |) -([28 | simplify considerable 

□E = 0, (35) 
(3) R ab = ^V Q EV h E. (36) 

We have pointed out that the rescaling ([31) simplifies 
the equations and allow us to write them in a more geo- 
metric form. This is the reason why this scaling is used 



in the case of U(l) cosmologies where the equations are 
locally the same but the norm rj never vanishes (see [TtJ 
[li ] and the review article [HI]). In our case the confor- 
mal scaling (|3"Tj) is singular at the axis. However, since 
the behavior of rj at the axis, as we will see in the next 
sections, is controlled a priori this singular scaling does 
not seems to introduce any extra difficulty in the equa- 
tions. We also remark that in all the numerical works 
mentioned above this conformal rescaling was not used, 
the equations are written in terms of the metric h ab de- 
fined by PD|) . 

Eqs. (l27 j) -([29 |) are purely geometric with respect to the 
metric h ab . To solve these equations we need to prescribe 
some gauge for the metric h ab . This will be done in the 
next two sections. 



B. 2+1 decomposition 

In order to formulate an initial value problem, we will 
perform an standard 2 + 1 decomposition of Eqs. {27J- 
(|29| . Note that this is completely analogous to the 3 + 1 
decomposition of Einstein equations, in fact all the for- 
mulas are formally identical because the dimension do 
not appears explicitly in them (sec, for example, [20| . 

a. 

Consider a foliation of spacelike, 2-dimensional slices 
S of the metric h ab . Let t be an associated time function 
and let n a be the unit normal vector orthogonal to S 
with respect to the metric h ab . The intrinsic metric on 
S is denoted by q a b and is given by 



(37) 



(38) 



Kb = -n a n b + q ab . 

Define the density fj, by 

[i = 2^R ab n a n b + 

and the current J b by 

■h = -q c b n a ^R cai (39) 

where ^R = {3) R a bh ab denotes the trace of (3) i? ab . 
Then, using Eq. (|29|) we obtain 



1 



2 + \Dlo\ 2 ) 



J a = -— 9 (v'DaV + lo'D a lo) . 
2rj £ 



(40) 
(41) 



where Da is the connexion with respect to qAB- The 
prime denotes directional derivative with respect to n°, 
that is 



1 



7]' = n a V al = - (d t T] - (3 a Dav) 



(42) 



where a is the lapse and j3 is the shift vector of the 



foliation. The indices A, B, 



denotes two dimensional 
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indices on S. The constraints equations corresponding to 
P5| are given by 



{2) R-x ab xab + x 2 = p, 

D A xab - D B X = Jb, 



(43) 
(44) 



where is the Ricci scalar of qAB, Xab is the second 
fundamental form of S and \ its trace 



x = q ab xab- 



(45) 



We use the following sign convention for the definition of 

XAB 



Xab = -q c a V c n b = -~£ n q a b, 



(46) 



where £ denotes Lie derivative. The evolution equations 
are given by 

d t qAB = -ZaxAB + £pqAB, (47) 
dtXAB = £/3Xab - D A D B a + ar AB , (48) 

where 

tab = XXab + {2) Rab ~ {3) Rab ~ ZxacXb- (49) 

and 



Rab = -J-^(dAVdBri + OaujObuj). (50) 



(3) 



The evolution equations (|47j) — (j48|) and the constraint 
equations (|4"3")) - (|4"4"1) constitute a complete 2 + 1 decompo- 
sition of the 3-dimensional Einstein Eq. (|29[) . It remains 
to decompose the effective matter Eqs. (|2?|) — (|28 [) . This 
can easily be obtained using the decomposition formula 
(jAllj) for the wave operator □ and the definition of the 
metric q a b given by (|3T)l . The result is the following 



(51) 



-oj" + A„cu + Daw— — — + w'x = % {D a wD a t 1 - u'jf) , 
a rj z 

(52) 

where instead of (|27|) we have use (1321) . and A q is the 
Laplacian with respect to qAB, i-e. A q = D a Da- 

Finally, we mention that the line clement of the metric 
h a b takes the standard form 

h = -a 2 dt 2 + q A B(dx A + (3 A dt){dx B + (3 B dt). (53) 



C. Gauge 

In this section we describe the maximal-isothermal 
gauge. In particular we review the mass formula for this 



gauge (see [lpj for details). For the lapse, we impose the 
maximal condition on the 2-surfaces 



x = o. 



(54) 



Note that we are not imposing that the surfaces are max- 
imal in the 3-dimensional picture as in [k| . The later 
condition is the one generally used Q [y]7but the dif- 
ference is only minor. In particular the mass formula is 
positive definite for both conditions as we will see. The 
one used here appears to be natural with respect to the 
rescaled metric h a \,. Eq. (|54[) implies the following well 
known equation for the lapse 



ot{x AB XAB 



where 



Mi 



{3) R ab n a n b 



— (v n 



(55) 



(56) 



The maximal gauge (|54p can be, of course, imposed in 
any dimensions and it is not related at all with axial 
symmetry. In contrast, the condition for the shift is pe- 
culiar for two space dimensions. The shift vector is fixed 
by the requirement that the intrinsic metric qAB has the 
following form 



qAB 



e 2u 5 AB , 



(57) 



where Sab is a fixed (i.e. d t dAB = 0) flat metric in two 
dimensions. Then, using (|54|1 . we obtain that the trace 
free part of (|4T1) is given by 



2a X 



AB 



(£ 9 /3) 



AB, 



(58) 



where C q is the conformal Killing operator in two dimen- 
sions with respect to the metric qAB defined in Eq. (|A1[) . 
Equation (|58|) is an elliptic first order system of equations 
for (3 A . 

The elliptic Eqs. ({55)) and (J3SJ) determine lapse and 
shift for the metric h ab and hence fixes completely the 
gauge freedom in Eqs. P?) - (f2T)]) . This gauge has asso- 
ciate a natural cylindrical coordinate system (t, p, z) for 
which the metric Sab is given 



dp 2 + dz 2 



(59) 



and the axis of symmetry is given by p = 0. The slices S 
are the half planes . 

For the analysis of the equations it is of course impor- 
tant to write them explicitly as partial differential equa- 
tions in these coordinates. We will do this in the remain- 
der of this section. In general, due to the complexity of 
Einstein equations, the partial differential equations ob- 
tained in a particular gauge can be quite involved. In our 
case, however, the geometric nature of the gauge plus the 
symmetry reductions will provide a relative simple set of 
equations. 

We first present some useful definitions. We need to 
subtract from r\ the part that vanishes at the axis. We 
define the function a by 



n 



P 2 e° 



(60) 
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Due to the rescaling ([3"Tj) , the lapse a vanishes also at 
the axis, hence we define the normalized lapse a by 



a = pa. 



(61) 



From the regularity conditions presented in the next sec- 
tion we will see that it is useful to define the function q 
defined by 



u = logp + a + q. 



(62) 



We now proceed to write the equations. We begin 
with the evolution equations for a and u>. The evolution 
equation for a is given by (|5Tj) . Using the definition (f^Uj) 
and the conformal rescaling expression for the Laplacian 
161) we obtain 

3 2 V + ( 3) A«t + d A a d -^- - 2e 2 "(logp)" + 2^ = 



o 

(-e 2u (Lu') 2 + \dLu\ 2 )p- 4 e- 2a . (63) 
In the same way, from ([52")) we get 



ap 



3 2U Lj" + (3) Auj + d A UJ 



d A c 



a 



- {-e 2u Jri + d A ojd A r,) . (64) 

7] 

Where d A denotes partial derivatives with respect to p 
and z and all the indices are moved with respect to the 
flat metric Sab- In these equations the lapse a and the 
shift (3 A appear trough the prime operator defined in 

621). 

The momentum constraint <|44j) is given by 

(65) 



9 b xab = Ja, 
where J A — e 2u J A , that is we have 

e 2u 

J A = - (r]'d A r] + u)'d A ui) . 
2rj z 



(66) 



To obtain (|65p we have used the conformal rescaling of 
the divergence in 2-dimensions given by (|A8|) . The in- 
dices in Eq. ([6"5")) and in the rest of the article, are moved 
with the flat metric S A b ■ To avoid confusion, it is useful 
to introduce the following notation 



(3a = P a 5 ab , Xb=^XCb, X^ = S a ^6^ Xcd . 

(67) 

That is, we want to distinguish between, say, the covector 
(3 A = f3 A q A B used in the previous section and (3 A (see the 
discussion after Eq. (|A10p in the Appendix). 

The Hamiltonian constraint, Eq. (|43[) . is given by 



where 



•2u 



(3) A(7 + Aq = -1, 



\dco\ 2 



(68) 



(t/ 2 + co' 2 ) + \da\ 2 + ^- + 2e^x AB XAB 



(69) 



Let us consider the evolution equations for q A s and 
Xab ■ The evolution equation for the metric q A B reduces 
to 



2dtu = d A p A + 2(3 A d A u. 



(70) 



And the evolution equation for the second fundamental 
form xab is given by 



dtXAB = £/3Xab - Fab - aG AB - 2cxx A cXb ( 71 ) 

where F A b denotes the trace free part (with respect to 
6ab) oiD A D B a. Using Eq. (|A4[) ) we obtain 

Fab = d A d B a - -^Aa - 2d( A ad B )U + d c ad c u5 AB . 

(72) 

And G A b denotes the trace free part of ^>R A b, namely 

G AB = {z) Rab - \s A b {3] R CD 5 cd , (73) 

where ^R-ab is given by (|50|) . 

The equation for the lapse is given by 

Aa = a(e- 2u X AB XAB+e 2u p 1 ), (74) 

and for the shift we have 

{Cf3) AB = 2ae- 2u X AB , (75) 

where C is the flat conformal Killing operator defined by 
®. 

Using the identity (|A3|) , we can transform the the first 
order system of Eqs. ([75]) for the shift and for the momen- 
tum constraint (|65p in a pair of second-order uncoupled 
equations. For the shift, we take a divergence to Eq. ([75]) 
to obtain 

A/3 A = 2d B (ax AB e- 2u ). (76) 

For Eq. l[6"5")) we define the vector v A by 

Xab = C(v)ab, (77) 

and hence Eq. ([75)1 transform to 

Av A = J a- (78) 

The total ADM mass of the spacetime can be calcu- 
lated as a volume integral on the half plane K 2 ^ of the 
positive definite effective energy density (see [13) 



1 f 

m = — I e pdpdz. 

16 J R 2 + 



(79) 



Finally, we mention that for the twist free case (u = 0) 
the four dimensional spacetime metric has a simple 
expression in these coordinates, namely 



g = — T e-°dt'+ 
P 2 



((dp + (3 p dt) 2 + (dz + p z dt) 2 ) + p 2 e a d(f> 2 . (80) 
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D. Boundary conditions and axial regularity 

The boundary conditions at the axis in axial symme- 
try have been extensible analyzed in the literature 
[3], 0, Hi- They involve parity conditions in the p de- 
pendence of the different fields. That it, the relevant 
functions are either even or odd functions of p. In order 
to use these results in our setting, it is useful to write the 
relations of the quantities with respect to the rescaled 
metric h a b and the canonical metric h a b, since all the 
above mentioned articles work with the metric h a b- 

Using relation (j3Tj) we obtain for the 2-dimensional 
metric 



QAB = VQAB, 

and for the second fundamental form 

XAB = y/rj XAB + 7T QAB 
V 2 V 



(81) 



(82) 



where quantities with a tilde are written with respect to 
the metric h a b- We also have 



a = ^a, I3 A = P A 



(83) 



Using these relations and the results mentioned above it 
is straightforward to obtain the following behavior of the 
relevant variables 



?7, u>, a, u, q, cr, xppi fl z are even functions of p, (84) 



and 



Xpz j P P are odd functions of p. 



(85) 



Note that odd functions vanishes at the axis and the p 
derivative of even functions vanishes at the axis. It fol- 
lows that one can impose homogeneous Dirichlet bound- 
ary conditions at the axis for odd functions and homoge- 
neous Neumann boundary conditions for even functions. 
In addition, we have that the function q defined by 
should vanished at the axis 



0. 



(86) 



Since q is an even functions, from (|86[) we deduce that 
q = 0(p 2 ) near the axis. Finally, there is an important 
regularity condition which comes from the axial regular- 
ity of the 3-dimensional extrinsic curvature. Let us define 
the following quantity 



1 

w = — 
P 



VXpp 



(87) 



Then it follows that 



w = 0(p), 



near the axis. This is the equivalent of the regularity 
condition given in equation (50) in [14 1 adapted to our 
conformally rescaled metric. See also [7 1, 



The fall off conditions at infinity are the standard 
asymptotically flat ones. In particular we have 



lim a — 1, 

r — >oo 



(89) 



and 



a, (3 A = Oir- 1 ), X AB=0(r- 2 ), (90) 



as r — > co. 



IV. LINEARIZED EQUATIONS 

In this section we make a linear expansion around 
Minkowski of the Einstein equations in the maximal- 
isothermal gauge described in the previous section. Note 
that for Minkowski we have 



V = P > 



(91) 



and hence, due to the rescaling (|3"Tj) , the background met- 
ric h a b, given in coordinates by (|53p . is non-flat 

h = p 2 (-dt 2 + dp 2 + dz 2 ) . (92) 

The other background quantities are given by 

uj = 0, a = p, (i A = 0, xab = 0, (93) 

and 

u = lnp, g = 0, (7 = 0. (94) 

The Hamiltonian constraint and the equation for the 
lapse are non-trivial for the metric (|92[) . namely 



Aa = 0, An = — . 

P 2 



(95) 



Let us proceed with the linearization. For simplicity 
we will consider only the case ui = 0. The first step is to 
compute the lapse function. The right hand side of Eq. 
([74")) is second-order, then, using the boundary condition 
we obtain 



a = l. 



(96) 



That is, the maximal condition for the lapse is trivial at 
the linearized level. On the contrary, as we will see, the 
equation for the shift plays a crucial role. 

The next step is to compute the linearization of the 
wave Eq. (|63j) for a, we obtain 



p+ (3) Ao- = 0, 



(97) 



where dot means partial derivative with respect to t and 
we have defined 



p = d t a 



2pP 



(98) 
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In order to close the system we need an equation for 
(3 P . Using equation (f6"5)) and (f6"6")) for the momentum we 
obtain 



with 



d A XAB = J A 



J A = ~pd A p 7 



(99) 



(100) 



We define the vector field v A by Eq. ([77)) and then by 
Eq. ([75]) we obtain 

Av A = -pd AP . (101) 
From poop we deduce J z = and hence we get 

Av z = 0. (102) 
By the fall off condition ([9"0"[) , we obtain 

v z = 0. (103) 
In the following, to simplify the notation we set 

v = v p . (104) 

Eq. (fTUTj) reads 

Av = -p. (105) 
Using (|77p we also obtain 

X PP = d p v, x P z = d z v. (106) 
For the shift we have the equation 



~AB 

(Cf3) AB = 2^ — . 

P 



(107) 



Taking a divergence to this equation (or linearizing (|76|) 1 
we obtain 



3 A na IX AB 



A0 A = 28 E 



(108) 



Note that in (|108|) we get an equation for j3 p decoupled 
from j3 z . Using Eq. (|100p , ([M)l from this equation we 
get 



A(3 p = - 2 -( P +^ 
P\ P 



(109) 



Eq. (T1TJ9]) . together with fTUS]) and (JHTJ form a complete 
system for the variables v, a and /3 P . Alternative, using 
Eq. (|107p and (|99|) we can eliminate xab and hence also 
v. We get the following equation for [3 A 



d B (p(CP) 



AB\ 



-2pd A p. 



(110) 



Eq. dTTUJ), together with fTUS]) and (JST]) form a complete 
system for the variables a, (3 P and f3 z . 



There is however an important difficulty. The lin- 
earization of the regularity condition (|87[) - (f55|) is given 
by 

1 / d v\ 
w = — [p+— ), w = 0(p). (Ill) 
P\PJ 



Were we have used Eq. (|106p . From the set of equa- 
tions presented above, it is difficult to ensure that this 
condition will be satisfied. To enforce this condition we 
will write the equations in terms of different variables. In 
order to do that, we need first to compute the remain- 
ing equations, namely the Hamiltonian constraint and 
the evolution equation for the metric and second funda- 
mental form. Since e defined in <|69j) is second-order, the 
Hamiltonian constraint (|68[) is given by 

Aq=-^Aa. (112) 
The evolution equation for q is obtained from (|70p 



2 p 



(113) 



The evolution equation for xab is obtained linearizing 

Xab = 2d (A qd B) p - S AB d p q. (114) 

We can also write the evolution Eqs. (|114p in components 

X PP = d p q, X P z = d z q. (115) 

Using Eqs. (|106p we deduce the important relation 

v = q, (116) 

which only holds in the twist free case. This relation 
simplify the equations considerable. From Eq. (I114p we 
also deduce 



d xab = Aqd A p. 



(117) 



With these equations we can compute the time deriva- 
tive of w 



If. d p q 



(118) 



q = pw + p d p [!—) . (119) 



and hence the evolution equation for q is given by 

'(3» 



As a consequence, q satisfies the following wave equation 
q = Aq-^+pd p \{jY (120) 

We also have 

9 / f>„n\ 

(121) 



Ap p = 2 -(Aq-^ 
P \ P 
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Eqs. (|120p and (|121[) form a complete systems for the 
variables q and [3 P . A similar choice of variable was used 
in by 3 However, in our particular case (i.e. linear 

equation without twist) it is possible a further simplifica- 
tion, namely to use Eq. (|116[) and hence replace q by v in 
these equations and then integrate in time. In this way 
we obtain our main Eqs. ((T|) and @. The advantage of 
using v as a variable is that the mass integral has a sim- 
ple expression in terms of v given by (|19p . This formula 
for the mass is obtained expanding up to second-order 
the energy density (|6"9"|). using equations (|106[) to replace 
Xab and Eq. (|105[) to replace 77' ' jr\. We discuss this in 
more detail in the next section. 

The boundary conditions and for (3 P arise from 
the axial regularity condition (|85|) and the asymptotically 
flat fall off 150)) . Conditions ([5]) arise from the axial reg- 
ularity conditions for q given by (|84|) and (|86|) . The main 
advantage of Eqs. (JTJ and ((2|) is that they have build in 
the regularity condition (jllip as we will see in the next 
section. Let us mention that the non-trivial regularity 
condition (jllip is written in terms of v as follows 

w = ^(av-^-Y w = 0{p). (122) 

The component (3 Z , which does not appears in Eqs. (TTJ 
and ([2]), can be calculated using 

d p (3?-d z p z =2^, (123) 

9,/3 p + 9p/3 z =2^. (124) 
P 

Or, alternative, using Eq. (fl6|) which is is obtained tak- 
ing a derivative to Eqs. (|123j) - (1124l) . Finally, the four 
dimensional perturbation (|10J) is obtained using the line 
element and the background values flU]), (JM]), (jM]) . 



V. PROPERTIES OF THE LINEAR EQUATIONS 

In this section we analyze some properties of our linear 
equations (H|) and ©. We begin with the symmetries of 
these equations. The first symmetry is given by trans- 
lation in z. This is to be expected since the gauge fixes 
the axis (and hence there is no translation freedom in 
p), but we still have the freedom to chose the origin in 
the z coordinate. Then, if we have a solution v,(3 p ; the 
derivative d z v,d z (3 p is also a solution, since d z commute 
with all the differential operators because their coeffi- 
cients depend only on p. The same argument applies to 
time translations, which is the second symmetry of the 
equations. The third symmetry is scaling. Let s a posi- 
tive real number. For a given solution v(t, p, z) wc define 
the rescaled function as 

v B (i,fi,2)=v(-,£,-\ (125) 



where 

t p z , 

t=~, P=~, £=-. 126 

s s s 

And the same for (3 P . Then, v s define also a solution in 
terms of the rescaled coordinates. The mass rescales like 

to — * sm. (127) 

In order to understand the equations in a simpler situ- 
ation, let us first consider Eqs. (jTJ) and j2| ona bounded 
domain il which does not contain the axis. On f2 the 
coefficient of Eqs. JT]) and ^ are smooth. Eq. is 
an elliptic equation for f3 p if we consider v as a given 
function. Hence, in order to solve this equation we need 
to prescribe elliptic boundary conditions for f3 p on dfl. 
For example, Dirichlet or Neumann boundary conditions. 
Eq. (fT]) is a wave equation for v if we consider f3 p as a 
given function. To solve this wave equation we need to 
prescribe initial data for v and v at t — together with 
compatible boundary conditions for v at <9f2. For example 
Dirichlet, Neumann or Sommcrfcld boundary conditions 
for v at dil. The equations are of course coupled, so it 
is not obvious that the above procedure of fixing bound- 
ary conditions is correct since v and f3 p are not "given 
functions". However, it is possible to prove that this is 
procedure is in fact correct. Consider the following iter- 
ation scheme 

v n+1 - Av n+1 + = pdpl — , (128 

p V p J 

Aft +1 = l(Av n -^y (129) 

In this iteration, the equations are not coupled and hence 
the boundary conditions mentioned above (which are 
kept fixed) are correct. Following similar arguments to 
the one presented in [22j (see also [23[) it is not diffi- 
cult to see that this iteration converges for some small 
time interval. And hence we get well-posedness for the 
linear system (|T|) and ^ under these boundary condi- 
tions on the domain f2. The reason why the iteration 
(|128p and (|129p converges is the following. From Eq. 
(|129p . using standard elliptic, estimates we obtain that 
(3 P is equivalent (in number of derivatives) to v. Hence, 
the term containing j3 p in Eq. (|128p is equivalent to a 
first order derivative of v and then it is not in the prin- 
cipal part of the wave equation. This rough argument 
suggests that the combination of elliptic estimates and 
energy estimates for the wave equations will close and 
hence the iteration will converge. This is basically the 
argument presented in [z| and [23| . If the domain Q 
is not bounded, this argument will still work if we add 
appropriate fall-off conditions at infinity. However, the 
situation change drastically when f2 includes the axis. 
Let us analyze that case. 

Since the axis is a singular boundary for the equa- 
tions, we are not free to chose arbitrary boundary con- 
dition there. In fact j3 p and d p v should vanishes at the 
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axis, otherwise the equations become singular. If we use 
L'Hopital rule, we conclude that the term with (3 P in ((T|) 
contain in fact two derivatives with respect to p at the 
axis. That is, due to the L'Hopital limit, to divide by p 
is equivalent as to take a derivative with respect to p at 
the axis. But then, using Eq. j2]), we conclude that this 
term is equivalent to second derivatives of v and hence it 
is in the principal part of the wave equation. We can not 
conclude that the iteration scheme (| 128|) — (| 129[) converges 
if we include the axis in the domain. This is, roughly 
speaking, the main difficulty to prove the well-posedness 
of the linear system (TT)) and @ . It appears to be diffi- 
cult to identify the principal part of the system at the 
axis and to construct an appropriate iteration scheme. 

Let us discuss in detail the boundary conditions at the 
axis. We are interested in solutions v which vanish at 
the axis, this comes from the regularity condition (|86p . 
Moreover, we have seen in section UlI D[ the smoothness 
of the spacetime metric at the axis implies that the func- 
tions (3 P and v should satisfies the parity conditions ([54")) 
and ([55)1 . Let us see heuristically, how these conditions 
are automatically implied by the equations provided we 
impose the following standard boundary conditions. At 
the axis we impose 

l3 p \ P =o = 0. (130) 
For v we prescribe initial data 

v\ t= o = /, v\ t =o = g, (131) 

such that 

/| P =o=0, 5 Ip=o = 0. (132) 

Note that we are not imposing any condition on v at 
the axis for t > 0. We make a formal series expansion, 
namely let as assume that our solution is smooth at the 
axis and has the form 

oo oo 

« = £p n On(i,«), /3 p = 5> n &„(i,z). (133) 

n=0 n=0 

Substituting these expansions in Eqs. ([1]) and ([2]) we ob- 
tain the following recurrence relation for the coefficients 

a n = (n + 2)na n+2 + d 2 z a n + nb n+ i, (134) 

and 

(n + l)n6 n+ i + d 2 z b n - x = 2(n + 2)na n+2 + d 2 z a n . (135) 

These expressions are valid for all integer n, with the 
convention that the coefficients b n and a n vanished for 
n < 0. The first non-trivial n in Eq. (| 1 34[) is n = — 1, 
which gives the relation 

oi + b = 0. (136) 

The term n — is given by 

do = d 2 ao- (137) 



This is a wave equation in 1-dimension. From the bound- 
ary conditions (|130[) we obtain 

b a = 0. (138) 
Hence we deduce from ([1361) that 

ax = 0. (139) 

From the initial data conditions (|132p we have that 

a |t=o=0, d | t=0 = 0. (140) 

These provides trivial initial data for the wave Eq. (|137[) 
and hence we deduce 

ao = 0. (141) 

That is, we have deduced the behavior v = 0(p 2 ) only 
from the boundary conditions (|130|) and the condition on 
the initial data (|131[) . We want to prove now that Q138|) 
and (|139[) imply that all a n with n odd and all b n with n 
even are zero. We prove this by induction. Let us assume 
that for some n (with n > 1) we have that 

6„_l = 0, a n = 0. (142) 

Using Eq. (|135[) we deduce 

(n + l)6„+i = 2(n + 2)a„ +2 , (143) 

and from (|134p we have 

(n + 2)a„ +2 = -b n+ i. (144) 

And then we have a n +2 = &n+i = 0. Since (|142p is valid 
for n = 1 we have proved the desired result. That is, the 
solutions v and f3 p satisfy the parity conditions ([54"| and 
(|85| respectively. Using that v is an even function of p 
and that v — 0(p 2 ) it is straightforward to deduce that 
the regularity condition (|122[) holds for all times. 

We analyze the fall off behavior of the solution v. This 
behavior is completely determined by the initial data / 
and g. Let us assume that the initial data has compact 
support. In the case of the wave equation, the signal 
will propagate with finite speed and hence the solution 
will always have compact support for any finite time. 
In our case, however, the coupling with the elliptic Eq. 
@ produce a non-local behavior. Even if we start with 
compactly supported data, the function (3 P will instan- 
taneously spread to all space. Let us perform a formal 
expansion in r to see the typical behavior of v. We have 
that f3 p = 0{r~ r ) for all times, this is prescribed by the 
boundary conditions. In [1 01 ) it has been proved that this 
implies that j3 p / p — 0(r~ 2 ), and hence the terms con- 
taining (3 P in (D) is 0{r- 2 ). Then, at t = we obtain that 
v = 0{r~ 2 ). If we take time derivatives of the equations 
and repeat this argument, we get that all time derivatives 
of v are 0(r~ 2 ). Then, we conclude that the typical fall- 
off behavior for asymptotically flat solutions is given by 
©, in the sense that we can not expect a faster decay 
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in general. Instead of compactly supported data we can 
begin with initial data for v such that they are 0(r~ 2 ) 
at infinity. 

Let us discuss now the most important property of 
equations ([I]) and J2]) namely the mass conservation. As 
usual, the mass appears as a second-order quantity that 
can be calculated in terms of squares of first order quan- 
tities. The density (|69p up to this order is given 



(Av) 2 + \daf 



(145) 



The total mass is calculated by the integral (|79p . The 
mass integral is conserved for the full nonlinear equa- 
tions in this gauge (see [l(|) and hence it is conserved at 
the linearized level. It is however important to compute 
explicitly this conservation formula using only the linear 
Eqs. H]) and Note that the function ct appears in 
the mass and this function should be calculated from v 
using Eq. Q12[) . To compute the time derivative of m we 
need first to calculate the time derivative of ct. Using the 
evolution Eq. ([1]) only, we compute 



( 3 )A(-A« + 2^) 
P 



where we have defined 



-Av - 



1 



d p [p(A/3P-2L(v))}, 



L(v) = d J 



Av — 

P 



(146) 



(147) 



Then, using ([2|) and the time derivative of (|12|) we get 



< 3 >A(-Av + 2^ -ct) = 0. 
P 



(148) 



If we are solving in the whole half plane then, by 
the fall-off conditions, we deduce that the only possible 
solution of this equation is the trivial one, and hence 



ct = -At 



Bp 



(149) 



We have proved that Eqs. {T]) and ([2|) together with 
(fT2|) imply Eq. (|149|) . We can also formulate the system 
in a different way. We can take fT} and © and Eq. (|149[) , 
instead of (fl"2|) . as an evolution equation for ct. If we take 
the Laplacian ( 3 ) A to both sides of Eq. Q149P and use the 
identity (|146[) together with Eq. @ we obtain 



( 3 ) Act = -Av. 



(150) 



Hence, if we chose initial condition for ct such that 

(3) ACT| f=0 = -A*| t=0 , (151) 



Eq. (|150[) implies (fl"2"|) . This two different ways of calcu- 
lating ct correspond to a constrained system and a free 
system (using the terminology defined in 0]). The pre- 
vious calculation is nothing but the propagation of the 
Hamiltonian constraint at the linearized level. For the 



full Einstein equations, the difference of constrained and 
free evolution schemes involves different set of evolutions 
equations. In our linear system the evolution equations 
are the same (namely, (fTJ) and @), the difference is the 
way the function ct (and hence the mass) is calculated. 
These two ways are of course completely equivalent when 
the domain is the whole half plane R^_, however, as we 
will see, they are not equivalent for a bounded domain. 

Using Eqs. fT4"9"|) . ([T2|). ([!]) and ([2"]) we obtain the fol- 
lowing local conservation law for the density e defined by 



pe = d A e / 



(152) 



where 



,d A v . 



e A = 8— « + 2pad A a + Af3 p d A v - 4vd A (3 p . (153) 

The vector e A can be interpreted as the energy flow of 
the gravitational field. If we integrate Eq. (|152|) in R^_ 
we have that the boundary terms vanishes both at the 
axis (by the axial regularity) and at infinity (by the fall 
off conditions). Then we have 



0. 



(154) 



We can also integrate Eq. (|152[) on a bounded domain 
f2, namely we define the mass contained in f2 by 



nin = I epdpdz, 



and then we have 



e A n A , 



no. 



(155) 



(156) 



where n A is the unit normal of dil. The quantity e A n A 
measure how much energy is leaving or entering the do- 
main. The local conservation formula (1 1 52[) can be gen- 
eralized for the non-linear equations [24|. 

Using the conservation of the mass (|154|) we can prove 
uniqueness of solutions of the system. Let us say we have 
two different solutions with the same initial data. We 
take the difference between the two solutions. The dif- 
ference satisfies the same equation with zero initial data. 
In particular ct on the initial surface is zero. And hence 
the mass is zero. Since it is conserved the mass is zero 
for all times, which implies that the solution is zero. 

In the case of hyperbolic equations (the wave equa- 
tion for example) the conservation of the energy gives 
also local properties of the solution, namely finite speed 
propagation of signals. However this is not the case here; 
the elliptic equation implies a non-local behavior of the 
solution. 

The discussion above applies for the domain which 
is the relevant domain for the equations. However, in nu- 
merical computation we need to solve the equations on a 
finite grid and hence it is necessary to impose boundary 
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conditions on a bounded domain. A typical domain for 
the numerics is shown in Fig. [TJ As we mention in sec- 
tion [TTJ for our present purpose we only need to prescribe 
some boundary conditions compatible with asymptotic 
flatness. For example, homogeneous Dirichlet boundary 
conditions for v and (3 P . However, the mass formula rise 
an interesting point here. On a bounded domain, to cal- 
culate a we have two possibilities. First, we can deter- 
mine a as the unique solution of the elliptic equation 
(fT2")) with some boundary conditions. If we do so, then 
we again deduce Eq. (|148p . However, from this equation 
we can not deduce (|149[) . In effect, we have 

P 

& + Av - 2— = H, (157) 
P 

where H satisfies 

(3) Aff = 0. (158) 

We can not conclude that H is zero from this equa- 
tion, because H will have non-trivial boundary condition. 
Namely, let us assume the we prescribe some boundary 
condition for a. We can not control the boundary value 
of Ad, and hence we can not ensure that H vanishes at 
the boundary. In fact, the function H is fixed as the 
unique solution of (| 1 58[) with boundary values 

H\en = {a + Av - 2^)| an . (159) 
P 

Then, if we compute the time derivative of the density e 
we get 

pe = d A e A + pHAv. (160) 

That is, we do not get a conservation law, there is a 
volume term given by H . There seems to be no boundary 
conditions for a that can ensure H to vanishes. 

The other possibility is to compute a using the evo- 
lution equation (|149[) with initial condition (|151[) . From 
(|149|) , in the same way as we mentioned above we deduce 
(I150p . since in this deduction the boundary conditions 
play no role. Using the initial data condition (|15ip . from 
(I150p we deduce (|12j) . That is, we are in the same situa- 
tion as the whole domain. Hence, in this case we recover 
(|152l) . where e A is given by the same expression (I153|) . 
From this point of view, this evolution scheme appears 
to be better than the previous one. 

In this scheme, we are free to chose any elliptic bound- 
ary condition for p p and any boundary condition for v 
compatible with the wave equation. For a we do not 
have any freedom, and hence we can not prescribe the 
boundary value of this function. 

A natural choice of boundary conditions would be to 
force the boundary integral in (|156p to have a definite 
sign. These conditions would have the interpretation of 
radiative boundary conditions, in the sense that the en- 
ergy is leaving the domain. To prescribe such conditions 
seems not to be possible (at least for generic data) since 



L 















P 


R 



FIG. 1: The bounded domain Q for the numerical evolution 



we do not have any control on the term with a. However, 
we can do something intermediate. Namely, if we impose 
Sommcrfeld boundary condition for v 

v = -n A d A v, (161) 

and homogeneous Dirichlet conditions for P we have 
that the first term in (|153[) has negative sign, the third 
term is zero. For the second and fourth term we have no 
control a priori. But we can expect that the influence of 
these term is small at least for some class of initial data. 
If this is true, then we get 

m n < 0. (162) 

This is what we observe in our numerical simulations 
described in the next sections. 



VI. NUMERICAL IMPLEMENTATION 

In this section we want to study numerically the initial- 
boundary value problem (IBVP) for the Eqs. {]]) and 
@. In this problem the symmetry axis, p = 0, becomes 
a boundary of our domain. Notice then, that working 
with the variable v poses an inconvenient as regards the 
boundary condition at p — since, according to (|8j). this 
function satisfies both, homogeneous Dirichlet boundary 
condition and homogeneous Neumann boundary condi- 
tion. It is then convenient to rewrite the equations in 
terms of a new variable for which the smoothness prop- 
erties at the symmetry axis defines a unique, equivalent, 
boundary condition. We define v — v/p. This new vari- 
able vanishes linearly with p and the correct boundary 
condition is simply homogeneous Dirichlet at p = 0. 

The equation for v(p,z,t), with p e [0, R], z G [0,L], 
and t > 0, is 

v = AC + d p Q +d P {fj), (163) 
where (3 p (p,z,t) is determined by the elliptic equation 
Ap? = 2 (^Av + d p (~j) (164) 
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with homogeneous Dirichlet boundary conditions, 



P p (0,z,t) = f3 p (R,z,t) = 0, ze[0,I], 
P p (p, 0, t) = fi p {p, L, t) — 0, pe[0,i?]. 



(165) 



for all t £ [0, oo). 

The boundary condition for v at the symmetry axis is 



v(0,z,t) = 0, 



(166) 



while at the outer boundaries we study two possibilities, 
homogeneous Dirichlet, 



v(R, z, t) — v(p, 0, t) = v(p, L, t) = 0, 

or Sommcrfeld (outgoing waves) 

v(R, z, t) = —d p v(R, z, t), 
v(p,0,t) = v z (p,0,t), 
v(p,L,t) = -v z (p,L,t). 



(167) 



(168) 



The initial data are 



v(p,z,0) = vo(p,z), 
v(p,z,0) = v t(p,z), 



(169) 



where vq and vot are C°° functions with compact support 
in (0, R) x (0, L) so that the compatibility of the boundary 
and initial data is not an issue. 

The Eqs. |T63)) - p69| constitute the IBVP we approx- 
imate with our finite difference scheme. 

We want to emphasize here an important difference be- 
tween our numerical approach with the usual approaches 
in the area (see for example [5j). We solve the IBVP 
for (|163p as a second-order equation just it is written 
above, i.e. we do not reduce (|163j) to a first order sys- 
tem of equations. The treatment of evolution equations 
as second-order equations as opposed to first order sys- 
tems of equations has several advantages. For example, 
the number of dynamical fields, and then the number of 
equations, is not increased. This facilitate the treatment 
of the boundary conditions. There are also numerical 
accuracy advantages. In the context of general relativ- 
ity, this has been stressed in (25j. In particular the sim- 
plest proofs of well-posedness for general initial-boundary 
value problems for Einstein's equations have been found 
recently using second-order systems of equations , [13] • 

a. The Implementation. In our numerical experi- 
ments we always consider square domains, i.e., R = L. 
To define the numerical grid let N be a positive integer 
and h = L/N the space stepsize. We define our grid to 
be half a stepsize displaced from all the boundaries. We 
think of our grid as a uniformly distributed set of points 
each of which is at the center of one of the N 2 square cells 
covering the domain. The coordinates of the gridpoint at 
the site (i, k) are then 



Pi 

Zk 



M< -3/2), 
h{k- 3/2), 



* = 0,1,2,.. 
fc = 0,l,2,. 



. N 



(170) 
(171) 



The sites (i,k) with 2 < i,k < 7V+1 are within the 
domain, while the sites with i = 0, 1, N + 2, JV + 3 and 
k = 0, 1, N + 2, N + 3 are "ghost points" used to ease the 
implementation of the boundary conditions Q- Time is 
discretized as 



t n = nSt, n = 0, 1, 2, 3, 



(172) 



We use capital latin letters to denote the grid functions 
associated to the dynamical variables. Also, we use sub- 
indices to denote the space-site indices and a super-index 
to denote the time step. This is, 



VTfc corresponds to v(pi, Zk,t„), 
B" k corresponds to (3 p (pi, Zk,t n ), 



(173) 



Besides the uniform grid we introduce the extra grid- 
points placed at the physical boundary 

(p = L,z k ), (pi,z = 0), (p l ,z = L) 

and denote the values of v at these points as 

VE,i, VP , VP L (174) 

respectively. 

In our difference scheme we approximate space deriva- 
tives by the standard fourth order accurate centered dif- 
ference operators given by [28] 



D := D ( 



D 



(175) 



and add a sub-index p or z to indicate what coordinate 
the operator is acting on. For example d\ v(pi, Zk, t n ) is 
approximated by 



/./•' 
-V; 



i,k— 2 



i,k—l 



i6 v: 



i,k+l 



k+2 



12h 2 

At every time step we need to solve the elliptic Eq. 
(|164p which we approximate by 



(Di + Di)B? tk =2[(Dj + Di)V& + 

% . k — 2,3,. 



D, 



Pi 

.N + l 



(176) 



We solve this difference equation iteratively using the 
Gauss-Seidel iteration scheme, and stop the iteration 
when the difference between both sides in (|176|) is 
smaller, in maximum norm, than a given small tolerance 
e. We then extend the solution to the ghost points — so 
that the homogeneous boundary condition is satisfied — 
as follows 



u o,k 

r>n 

n N+2,k 
on 



l.fc) 



- B 3,ki 
nil 



nil 
~ n i,N+li 



rtn 

n N+3,k 
qn 

n i,N+3 



pn 

pn 
~ n N,k 

-B? <2 

pn 

i.N 
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We now describe how we approximate (| 1 63[) using 
fourth-order accurate difference approximations in space; 
to use second-order accurate approximations instead, we 
just need to change D and D 2 in what follows by Do and 
D + D- respectively. 

t = 0. We set 

V° k =v (pi,z k ), i,k = 2,3,...N + l (177) 

and extend the solution to vanish at all ghost points 
and boundary points since the initial data has compact 
support. Then we compute B® k by solving (| 176[) as ex- 
plained above. 

t = St (first step). We do, for i, k = 2, 3, . . .N + 1, 



V£ k = VP k + St v ot ( Pt ,z k ) + -(St) 2 ((Dp + Dl)V» k 

+D p (V? >k / Pi )+D p {Bl k / Pi )) (178) 



Now, if working with boundary condition (|166[) , (|167p , we 
define the solution at the boundary points to vanish 



We notice that the second derivative in time is approx- 
imated by D + D_ which is second-order accurate. The 
time step we use in all our runs is St = h/ 10. The ra- 
tio St/h = 0.1 satisfies the Courant condition and we see 
from our runs that the whole method turns out to be 
numerically stable. 

Besides the solution V" k and J3" fc , an essential quan- 
tity we want to compute is the mass m^(t), defined by 
(|145p and (|155p . during the whole evolution. To this end 
we need to compute a(t) on the physical domain at all 
times. Given the approximations of v and (3, we compute 
a(t) by integrating (|149p — rewritten in terms of V, as an 
ODE at each gridpoint. The initial data for these ODEs 
is computed by solving the elliptic Eq. (|15ip , also rewrit- 
ten in terms of v, only once at initial time with homoge- 
neous Dirichlet boundary conditions and using the same 
technique we use to compute (3. The first time step to 
integrate (|149p is carried out with explicit Euler method, 
and from there on with the two-step, second-order ac- 
curate, Leap-Frog method. We evaluate the integral in 
(|155p with the midpoint rule. 



vh, k = v£ = vIl = o, 



(179) 



while if working with boundary condition (| 1 661) . (| 1681) we 
evolve the boundary points by integrating the boundary 
condition using explicit Euler scheme. For example for 
the boundary p = L 



vL = v?k-stD"vl k , 



(180) 



where 



24h 



is a fourth-order accurate approximation of the normal 
first derivative at the border p = L. We now extend the 
solution to the ghost points as 



v ) k 


= -vi, k , 




vl k 


= -v 2 ) k 






= 2V2, fc - 




^W+3.fc 


= - 












= 2^ - 






= 2^ L - 




Vi,N+3 


= - 





Finally, we compute Bj k as explained above. 

At t = n St. With n — 2, 3, . . . we evolve the solution 
with the two step method 

+D p (V^- 1 / Pi ) + D^B^/Pi)) (181) 

for i, k = 2, 3, . . . JV + 1. Then we impose the boundary 
conditions exactly as done in the first step. Finally we 
compute B" k as explained above. 



VII. TESTS, RUNS AND NUMERICAL 
RESULTS 



The numerical calculations we carry out in this work 
pursue two main objectives. The first objective is to 
make plausible that the initial-boundary value problem 
for (|163p . (|164p is well-posed. If we were simulating an 
IBVP that is not well-posed, the expectation would be 
that almost any consistent numerical simulation of the 
problem would fail to pass convergence tests, numerical 
stability tests, or both. We show below that both kind 
of numerical tests are passed satisfactorily by our nu- 
merical approximation. The second objective is to study 
the behavior of the mass in these initial-boundary value 
problems. In particular, we will show that for fixed ini- 
tial data, the larger the domain used in our calculation 
is, the longer and better the mass approaches a constant 
value. 

We use in our runs two kinds of initial data which 
are smooth and strongly decaying outside a small region 
(Gaussian functions). The first is 



v Q (p,z) = exp( 
vot(p,z) = 0. 



(p- 1/2) 2 + (z-L/2f 
0.1 2 



(182) 



which decays very fast as (p, z) get away from (1/2, L/2) 
and so approximate very well a compact support data 
on the domains we use. The second is the same kind of 
function but for v ot instead of Vq. namely 



Vo(p, z) 
vot(p, z) 



0. 



50 exp 



(p- l/2f + (z~L/2f 
0.1 2 



(183) 



Linearity of the problem tells us that the runs with 
Vq 7^ 0, or VQ t 0, can be performed separately. A 
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solution with general initial data is the superposition of 
two solutions, one with each kind of data. 

b. Elliptic Solver Tolerance. We need to determine 
the value of e to use in our runs. To this end we perform 
runs for six different values of e with all other parameters 
fixed to typical values in our runs. In these tests runs we 
use initial data given by (|182p and Sommcrfeld bound- 
ary conditions (|166[) . (|168[) . We then analyze the different 
values of the mass obtained for the six solutions. By com- 
paring the variations of ran (t) with respect to the initial 
value of mo, we see from our runs show that the evolu- 
tion is not very sensitive to the tolerance e. Fig. [2] shows 
that the plot for the different computed masses superim- 
pose when plotted in the full mass scale. The detail in 
the figure shows convergence of toq as e — > 0. 





FIG. 2: Evolution of the mass for six different values of the 
tolerance e. The upper plot shows the six runs in full mass 
scale. In this scale the six curves look superimposed. The 
lower plot shows a detail of the initial "flat" region in ampli- 
fied scale. 

In table Q] we show the maximum absolute difference 
between the computed masses with respect to the most 
accurate one (corresponding to e = 10 -6 ) and the time 
of occurrence. 

Based on this test we choose to use e = 10~ 3 in our 
further runs, which gives more that necessary accuracy 
for our discussion (around 10 -6 relative error.) 



e 


Amn 


tmax 


5mn/mno 


10" 


1 


1.27 x 10" 


-3 


0.390 


1.19 x 10" 4 


10" 


2 


1.81 x 10" 


-4 


0.406 


1.69 x 10" 6 


10" 


3 


1.45 x 10" 


-5 


0.481 


1.36 x 10" 6 


10" 


4 


1.08 x 10" 


-6 


0.503 


1.01 x 10" 7 


10" 


5 


8.19 x 10" 


8 


0.509 


7.65 x 10" 9 



TABLE I: Values of Amn = maxt \mn e (t) — m ni0 -6(t)\ and 
time of occurrence for different values of the tolerance e. Ini- 
tial mass is rano = 10.7083. 



c. Convergence Tests. To study convergence of the 
numerical solution we perform two series of runs in a 
unitary square domain with the initial data (|182|) . In 
the first series we use homogeneous Dirichlet boundary 
conditions and in the second Sommerfeld boundary con- 
ditions. 

Each series consists of four runs. In the successive 
runs we use h = l/N with N = 50; 100; 200; 400. In all 
runs St = h/10. Thus, in the second, third and fourth 
runs both h and St are divided by 2 with respect to the 
previous run. Let us call V^ h '(t) the solution computed 
using mesh-size h. 

The first, simplest and indirect, convergence test is to 
plot the masses for each run as a function of time and 
check, graphically, whether they converge as the value of 
h diminishes. Figures [3] and 0] show that this is in fact 
the case. 

A second more strict convergence and accuracy test is 
as follows. We compute the L2 norm of the difference 
between two successive runs. A simple analysis shows 
that, when the method is convergent and the mesh and 
time-step sizes are small enough, the quotient 



Qh{tn) 



\\V(h) {tn) _ v (h/2) {tn)h 
\\V^){t n )-VWi){t n )\\ 



(184) 



approaches the value 2 P where p is the accuracy order 
of the method. Our method is fourth-order accurate in 
space and second-order in time. Therefore the expecta- 
tion is that we obtain values of Qh that are close to 4 at 
most times. 

To compute the L2-norms we use the midpoint rule to 
approximate the integration on the coarsest grid of the 
two solutions being subtracted. Notice that the coarse 
grid is not sub-grid of a the fine one, as they are displaced 
from the domain boundaries by different amounts. Then, 
to evaluate the finest solution on the coarse grid we need 
to interpolate this solution. To do this we use bilinear 
interpolation. 

The results of this analysis are shown in the table ILTl 
and IIII1 The test is passed satisfactorily. 

d. Stability Tests. Numerical stability means that 
the solution to the IBVP stays bounded during time evo- 
lution. Typical signs of instability are the appearance of 
artifacts in the plot of the solution as a consequence of 
evolution and in most cases, after a while, the complete 
break-down of the solution. If an instability has its root 
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N=50 - 
N=100 
N=200 
N=400 



h=1/50 - 
h-1/100 
h-1/200 - 
h=1/400 



N=50 - 
N=100 - 
N=200 
N=400 





FIG. 3: Mass as function of time for evolution with homo- 
geneous Dirichlet boundary conditions. In the upper plot, in 
full mass scale, the four curves look almost superimposed. In 
the lower plot a detail in expanded mass scale shows that the 
curves converge to a limit curve when h and St diminishes 



FIG. 4: Mass as function of time for evolution with Som- 
merfeld boundary conditions. In the upper plot, in full mass 
scale, the four curves look almost superimposed. In the lower 
plot a detail in expanded mass scale shows that the curves 
converge to a limit curve when h and St diminishes 



on the ill-posedness of the analytic problem underneath, 
the expectation is that some high frequency modes of the 
solution explode exponentially fast and are detected at 
very short times of the numerical evolution. For some 
more benign ill-posed problems (like weakly-hyperbolic 
problems) the growing of instabilities is only polynomial 
and it may take longer to detect them. 

We performed several series of runs using both kinds of 
boundary conditions p^ . p^T]) or (fll^) , (fll)gj) and both 
kinds of initial data (|182[) or (|183[) on different domains 
and during several time intervals. We studied the plots 
of the solutions in all cases and they always look smooth, 
agreement with the boundary conditions imposed and 
never showed any sort of strange artifact. Typical plots 
for v(p, z, t) are shown in Fig. [5] We have also studied the 
plots of (3{p, z, t) in these runs and no sign of instability 
showed. 

A second, physically meaningful, test for stability is 
provided by the study of the mass toq which in this 
problem is a sort of incomplete H 2 Sobolev norm of 
the solution. As explained in Sec. [V] the mass is con- 
served for the Cauchy problem in the whole space. On 



t 


||y(2) _y(8) 


\l 2 




\l 2 


Qh(t) 


0.05 


5.2744x10' 


-5 


1.3214x10" 


5 


3.9913 


0.10 


8.2053x10" 


-5 


2.0665x10" 


5 


3.9706 


0.15 


9.2690x10" 


-5 


2.3339x10" 


5 


3.9715 


0.20 


9.8158x10" 


-5 


2.4938x10" 


5 


3.9360 


0.25 


1.1260x10" 


-4 


2.8857x10" 


5 


3.9020 


0.30 


1.3325x10" 


-4 


3.4668x10" 


5 


3.8436 


0.35 


1.6185x10" 


-4 


4.4464x10" 


5 


3.6401 


0.40 


1.8421x10" 


-4 


4.9874x10" 


5 


3.6934 


0.45 


2.2492x10" 


-4 


6.9417x10" 


5 


3.2402 


0.50 


2.9719x10" 


-4 


6.9240x10" 


5 


4.2923 



TABLE II: Convergence and accuracy quotient for solutions 
with homogeneous Dirichlet boundary condition. On the 
coarsest grid h — 10 -2 . 



bounded domains this is no longer true, but we expect 
that it stays bounded when using homogeneous Dirich- 
let boundary conditions, and that it goes to zero when 
using Sommerfeld boundary conditions. We analyze the 
behavior of the mass below. 



LS 



t 


||y(2) _ y(3) 


\l% 


llv-( 3 ) - v (4) 


\l 2 


Qh(t) 


0.04 


6.6922x10" 


5 


1.6703x10" 


5 


4.0067 


0.08 


6.2007x10" 


5 


1.5756x10" 


5 


3.9355 


0.12 


9.1495x10" 


5 


2.2958x10" 


5 


3.9854 


0.16 


9.2694x10" 


5 


2.3357x10" 


5 


3.9686 


0.20 


9.8226x10" 


5 


2.4911x10" 


5 


3.9431 


0.24 


1.0935x10" 


4 


2.7940x10" 


5 


3.9138 


0.28 


1.2407x10" 


4 


3.2083x10" 


5 


3.8671 


0.32 


1.4427x10" 


4 


3.8153x10" 


5 


3.7813 


0.36 


1.6740x10" 


4 


4.6440x10" 


5 


3.6047 


0.40 


1.8389x10" 


4 


5.0471x10" 


5 


3.6434 



TABLE III: Convergence and accuracy quotient for solutions 
with Sommerfeld boundary condition. On the coarsest grid 
h=10~ 2 . 



e. Behavior of the Mass. As explained before the 
mass, defined by Q145P and (|155p . is a conserved quan- 
tity when the Cauchy problem is considered in the whole 
space (i.e., SI is K+). In our numerical tests we solve 
the initial boundary value problem on compact domains 
where no known boundary conditions imply mass con- 
servation. However, an interesting study for the mass 
evolution can be done as follows. We solve the IB VP on 
domains of different size but use, in all runs, the same ini- 
tial data, at the same distance from the symmetry axis. 
The initial data are chosen to decay exponentially fast 
outside a region which is small compared to the smallest 
of the domains we use. Clearly, the expectation is that 
the larger the domain is the closest to constant the mass 
stays during evolution. 

We do series of runs for homogeneous Dirichlet bound- 
ary conditions and for Sommerfeld (outgoing waves) 
boundary conditions. The plots for the Dirichlet case are 
shown in Fig. [(51 Observe that the plot is not on full mass 
scale. The three curves show an almost constant initial 
region and then variations of small relative amplitude. 
After an initial peak immediately after the constant re- 
gion the amplitude of the variations is, roughly speaking, 
2% for the 1.28 x 1.28 domain, 1% for the 2.56 x 2.56 do- 
main and 0.6% for the 5.12 x 5.12 domain. The amplitude 
clearly diminishes when the domain increases size. For 
the case of Sommerfeld boundary conditions, the plots of 
the mass evolution can be seen in Fig. [7l This series of 
three runs is totally analogous to the previous case. The 
only change is the boundary condition used. As can be 
inferred from the plot in full mass scale, the energy leaks 
though the boundary as expected. 




FIG. 5: Plots of the solution v(t). Both plots of solutions 
computed on a grid with 128 x 128 gridpoints and initial data 
given by (|182|) . Upper plot is the solution with homogeneous 
Dirichlet boundary conditions at time t — 3.0, while lower 
plot is the solution with Sommerfeld boundary conditions at 
time t = 1.25. 
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VIII. FINAL COMMENTS 

In this article we have deduced the linear system 
and ^ and we have analyzed some of its properties. 
Among them, the most relevant are the mass conserva- 
tion and the numerical stability. The main open problem 



FIG. 6: Evolution of the mass mn as a function of time for 
three solutions with homogeneous Dirichlet boundary condi- 
tions and the same initial data but on domains of different 
size. In the upper right corner the each curve is associated to 
the corresponding domain. 
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1 .28 x 1 .28 

2.56 x 2.56 

5.12 x 5.12 



FIG. 7: Evolution of the mass mo as a function of time 
for three solutions with Sommerfeld boundary condition, the 
same initial data but on domains of different size. In the 
upper right corner the each curve is associated to the cor- 
responding domain. The lower plot shows in amplified scale 
that the "flat" region presents very small variations of around 
0.03%. 



kind of approach to this problem; although, of course, 
always resticted to axial symmetry. The ultimate and 
difficult goal is to say something, in this gauge, about 
the non-linear stability of a black hole in axial symme- 
try. 

The second road is the study axially symmetric per- 
turbation but with a black hole as background solution. 
Linear stability of the Kerr black hole is a relevant open 
problem which is currently intensively studied (see the 
review articles [3(|, [3l[ and references therein). The 
expectation is that the mass formula can help to prove 
linear stability under axially symmetric perturbation of 
the Kerr black hole. 



Acknowledgments 

S. D. thanks Piotr Chrusciel and Helmut Friedrich for 
useful discussions. These discussions took place at the 
Mathematisches Forschungsinstitut Oberwolfach during 
the workshop "Mathematical Aspects of General Relativ- 
ity", October 11th - October 17th, 2009. S. D. thanks 
Andres Acena for useful discussions that took place at 
the Max-Planck-Institut fur Gravitationsphysik (Albcrt- 
Einstein-Institut) during the conference "Space, Time 
and Beyond", October 19th - October 21th, 2009. S. 
D. thanks the organizers of these events for the invi- 
tation and the hospitality and support of the Mathe- 
matisches Forschungsinstitut Oberwolfach and the Max- 
Planck-Institut fur Gravitationsphysik (Albert-Einstein- 
Institut). 

S. D. is supported by CONICET (Argentina). This 
work was supported in part by grant PIP 6354/05 of 
CONICET (Argentina), grant 05/B415 Secyt-UNC (Ar- 
gentina) and the Partner Group grant of the Max Planck 
Institute for Gravitational Physics, Albert-Einstein- 
Institute (Germany). 



is to prove that this system is well-posed. Remarkable 
enought, it seems to be not much literature on this class 
of linear systems which are singular at the axis. 

Once the well-posedness problem is solved, we believe 
that the future research on the subject can be divided 
in two paralel but complementary roads. The first one 
is to extend the well-posedness from the linear system 
to the full Einstein equations in the maximal-isothermal 
gauge. The non-linear lower order terms introduce extra 
difficulties (see [y]). There are many possible evolutions 
schemes (see the discussion in Q). It is very likely that 
few of them (or may be only one) are well-posed. If this 
is the case, the resolution of the well-posedness question 
will lead us to select (or even discover) the correct evolu- 
tion scheme. After the local problem is solved, the next 
step is to use the global conservation of the mass to con- 
trol the full non-linear evolution in this gauge. A natural 
first example would be to recover the non-linear stability 
of Minkowski [29j in this gauge. The expectation is that 
the mass formula will provide a simpler (and different) 



APPENDIX A: USEFUL FORMULAS 

We collect in this appendix some useful formula that 
are used in the main part of this article. The confor- 
mal Killing operator in 2-dimensions with respect to the 
metric q AB is defined by 

(C q P) AB = D A f3 B + D B p A - q AB D c (3 c . (Al) 

For the particular case of a flat metric S AB this definition 
reduce to 

(Cf3) AB = d A l3 B + d B (3 A - S AB d c (3 c . (A2) 

For this operator we have the following identity often 
used in the article 

d B {£(3) AB = A(3 A . (A3) 

The Christofcll symbols of the metric q AB defined by 
(f57|) are given by 

Y C AB = 5%d A u + 6^d B u - d c u5 ABl (A4) 
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and the Ricci tensor is given by 

WR A b = -Au5 ab , ^R=-2e 



'At 



(A5) 



Under the conformal rescaling (|57|) the diferential op- 
erators relevant in this article transform as follows 



C q {l3) AB = e 2u C0) 



ABi 
4u o «AB 



where we have defined 

(3 A = e 2u p A X 



A B _ —AuSjAB 
- e X 



(A6) 
(A7) 
(A8) 

(A9) 



We follow the convention that the indices for hat quan- 
tities are moved with the flat metric Sab and inidices of 



non-hat quantities with the metric q AB . Then, we have 

XAB=XAB, P A = P A . (A10) 

That is why we suppress the hat notation for the tensors 
Xab and [3 A in the main part of this article. 

Take an arbitray spacelike foliation on Af,h a f,. The 
2 + 1 decomposition of the wave operator is given by 

D A a 

nf = -f" + A g f + D A f— + f X , (All) 

a 

where have made use of the following useful formulas 

n-V n A = ^ ) n«V a n t = ^. (A12) 
a a 
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